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The mode-sum method provides a practical means for calculating the self force acting on a small 
particle orbiting a larger black hole. In this method, one first computes the spherical-harmonic 
Z-mode contributions Fj 1 of the "full force" field F M , evaluated at the particle's location, and then 
sums over I subject to a certain regularization procedure. In the frequency-domain variant of this 
scheme the quantities Fj* are obtained by fully decomposing the particle's self field into Fourier- 
harmonic modes Imco, calculating the contribution of each such mode to Ff 1 , and then summing over 

00 ■ uj and m for given I. This procedure has the advantage that one only encounters ordinary differential 

equations. However, for eccentric orbits, the sum over ui is found to converge badly at the particle's 
location. This problem (reminiscent of the familiar Gibbs phenomenon of Fourier analysis) results 

£Nj ■ from the discontinuity of the time-domain Fj* field at the particle's worldline. Here we propose a 

! * ' simple and practical method to resolve this problem. The method utilizes the homogeneous modes 

Imui of the self field to construct Ft (rather than the inhomogeneous modes, as in the standard 
method), which guarantees an exponentially-fast convergence to the correct value of Fj 1 , even at the 
particle's location. We illustrate the application of the method with the example of the monopole 
scalar-field perturbation from a scalar charge in an eccentric orbit around a Schwarzschild black hole. 
Our method, however, should be applicable to a wider range of problems, including the calculation 
of the gravitational self-force using either Teukolsky's formalism, or a direct integration of the metric 
' perturbation equations. 

^jrj; I. INTRODUCTION 

! The problem of calculating the gravitational self-force [H, 0, Q actingon a pointlike test particle as it moves in orbit 
around a black hole is attracting considerable attention in recent years [^Ja- Within this context, various authors have 
' been studying also the analogous problem of the scalar-field self-force Q acting on a particle endowed with a scalar 
charge, which proved to be a useful toy model. The electromagnetic (EM) self-force, acting on an electric point charge, 
was also studied by various authors following the seminal work by DeWitt and Brehme 7]. A practical algorithm, 
commonly used for computing the self force in all three cases (gravitational, EM and scalar), is the mode-sum method 
This method requires as input the multipole modes of the full (retarded) perturbation fields, along with 
their derivatives, evaluated at the particle's location. These multipoles can be calculated using either frequency- 
domain methods (as in, e.g., [Ill, [H, El ) or time-domain numerical evolution (as in, e.g., 0, E EB, [I3| ) • I n both 
computational approaches one sets off by writing down the appropriate set of perturbation equations, modeling the 
source term associated with the point particle (namely, the energy-momentum, the electric four-current, or the scalar 
charge) as a delta-function distribution. In the frequency-domain approach one then decomposes the inhomogeneous 
perturbation equations into Fourier-harmonic modes ( a lmu> modes") and proceeds by solving the resulting ordinary 
differential equations (ODEs) with suitable boundary conditions at spatial infinity and at the event horizon. In the 
alternative, time-domain approach, one refrains from decomposing the field into frequency modes, and instead tackles 
the partial differential equations for each l,m directly using time evolution. 

Each of the above two approaches has its strengths and weaknesses. The time-domain approach has the advantage 
that one only deals with a single field for each Im, whereas in the frequency-domain approach one has to sum over the 
various u modes. On the other hand, the latter approach has the obvious advantage that one only faces ODEs. Despite 
the fact that time-domain methods are winning growing popularity in recent years, frequency-domain calculations 
remain an appealing option for some range of orbital parameters [l8j . Also, it turns out that the non-radiative 
multipoles of the gravitational perturbation in Schwarzschild (i.e., the modes 1 — 0,1) are difficult to analyze in the 
time domain, due to instabilities (l5l . [l7l | , and one resorts to a frequency-domain calculation in this case (l9j . Working 
in the frequency domain, however, brings about a technical issue which, to the best of our knowledge, has not been 
addressed so far in the current context. 

To illustrate the problem, is it instructive to refer to the simple case of (minimally-coupled and massless) scalar- 
field perturbations from a pointlike scalar charge orbiting a Schwarzschild black hole. In this case, the scalar field 
&(t,r,8,ip) can be decomposed in spherical harmonics Yi m (9,ip), yielding the multipolar mode functions <fii m (t,r). 
Here t, r, 8, <p are the standard Schwarzschild coordinates. Let us denote the r value of the particle's location at time 
t by r p (t). With suitable boundary/initial conditions, a unique solution is obtained for tfii m (for each l,m), which is 
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continuous along r = r p (t). However, the derivatives (j)i m ,r and 4>i m ,t will generally suffer a finite jump discontinuity 
across r = r p (t), which reflects the presence of a source "shell" resulting from decomposing the point charge in 
spherical harmonics. In particular, if the orbit is eccentric, the derivatives of <fii m will generally be discontinuous 
functions of t at a fixed value of r along the orbit. 

Now imagine trying to reconstruct 4>i m (t, r) (for some fixed r along the orbit) as a sum over its Fourier components: 



Since, for an eccentric orbit, <pim(t^) is only a C° function of t at the particle's worldline, it follows from standard 
Fourier theory 20] that the Fourier sum in Eq. |T]) will only converge there like ~ lo^ 1 . The actual situation is even 
worse, however, because for self-force calculations we need not only <fii m (t,r) but also its derivatives. For instance, to 
calculate the r component of the self force we need to evaluate (one-sided limits of) </>i m ,r(t, r) as r approaches r p (t). 
Suppose again that we want to reconstruct this quantity from its Fourier components (namely the functions Ri m u,r)- 
Since 4>im,r{tj r ) is a discontinuous function of t, we will inevitably face here the well known "Gibbs phenomenon" 
pH | . Namely, the Fourier sum will fail to converge to the right value at r — » r p (t) . (The problematic behavior of the 
Fourier sum is simply a consequence of our attempt to construct a discontinuous function — or a non-smooth function 
in the case of the field 4>i m itself — from a sum over smooth harmonics.) 

From a practical point of view this would mean that (i) at the coincidence limit r — > r p (t) the sum over to modes 
would fail to yield the correct one-sided values of (f>i m)r (t,r), however many lo modes are included in the sum; and 
(ii) if we reconstruct 4>i m ,r at a point r = tq off the worldline, then the Fourier series should indeed converge; Alas, 
the number of lo modes required for achieving a prescribed precision would grow unboundedly as tq approaches r p (t) , 
making it extremely difficult to evaluate <fiim,r at the coincidence limit. 

This technical difficulty is rather generic, and will show also in calculations of the local EM and gravitational fields. 
Consider, as a second example, the gravitational perturbation from a point mass moving in an eccentric orbit in 
Schwarzschild: In suitable gauges (like the Lorenz gauge, often applied in self-force calculations) the multipole Im- 
modes of the physical metric perturbation [26j ] are again C° functions of r and t along the orbit, and their derivatives 
are generally discontinuous there. Attempting to construct them naively from a sum over frequency modes would 
encounter the same difficulty as in the scalar case: A poor convergence of the (Zm-modes of the) metric perturbations, 
and lack of convergence for their derivatives. 

For orbits in Kerr spacetime the situation is basically similar though more subtle. The Kerr variant of the mode- 
sum method [lfj requires, just as in its Schwarzschild counterpart, the spherical harmonic modes <fii m (r,t) of the 
perturbation field (as well as their derivatives) as input. [13] For given Im, the spherical-harmonic decomposition of 
the point charge will again result in a <5-function-type source term distributed over a shell, which in turn renders the 
derivatives of (f>i m (r, t) discontinuous. Therefore an attempt to construct (f>i m (r, t) (and, more crucially, its derivatives) 
through a naive summation over its Fourier modes will lead to the same difficulties as in the Schwarzschild case. 

The problem discussed here takes an even more extreme form when considering EM or gravitational perturbations 
via the Teukolsky formalism: Here, the Im modes of the perturbation fields (now the Newman-Penrose fields tpo, Lp 2 or 
^o, ^4) are not even continuous at the particle's orbit — a consequence of the fact that the source term for Teukolsky's 
equation involves derivatives of the electric four-current or the energy-momentum tensor associated with the particle (a 
single derivative in the EM case; a second derivative in the gravitational case). 28] Again, a naive attempt to construct 
these multipoles as a sum over their lo modes will be hampered by the Gibbs phenomenon, and the associated lack of 
convergence. 

In this article we propose a way around the above problem, which is both elegant and extremely simple. In our 
method we use the homogeneous radial functions Rimui(f) (extended all the way through to the particle's worldline), 
instead of the actual inhomogeneous functions. The Fourier sum of these homogeneous radial functions is found to 
converge exponentially- fast, and to yield the correct values of the perturbation multipoles (and their derivatives) 
along the particle's worldline. We shall focus in this paper on the scalar-field case. We justify our new method using 
simple analytical arguments, and then demonstrate the validity of the method (and the exponential convergence) 
with an explicit numerical calculation in the case / = 0. The same method should be applicable, however, for any 
of the other problems mentioned above: EM and gravitational perturbations using Teukolsky's equation (or Sasaki- 
Nakamura's equation), as well as metric perturbations in the Lorenz gauge. A forthcoming paper [l9j will report on 
the computation of the local monopole and dipole modes of the Lorenz-gauge perturbation (for eccentric orbits in 
Schwarzschild), facilitated by the new method suggested here. 

This paper is structured as follows. In Sec. II we set up the physical scenario — a pointlike scalar charge in a bound 
orbit around a Schwarzschild black hole — and review the formalism commonly used in this case to construct the 
scalar-field multipoles and the scalar self-force. Section III demonstrates how the naive sum over frequencies leads 
to the Gibbs phenomenon and to the associated problematic convergence at the particle's location. Then in Sec. 
IV we present our new method of extended homogeneous solutions, and show how it cures the problematic behavior 




(1) 



3 



of the Fourier sum. We provide the theoretical justification to this method, as well as numerical verification in the 
monopole (/ = 0) case. In Sec. V we highlight the advantages of the new method and discuss foreseeable applications. 
Appendices lAUCl give details of the methods used for our numerical illustrations, and App. [Pi contains some technical 
details relating to the formal justification of our new method. 

Throughout this work we use standard geometrized units (with c = G = 1) and metric signature ( — h++). 



II. PRELIMINARIES 



A. Physical setup and scalar-field equation 

Consider a pointlike particle which moves on an eccentric, bound geodesic orbit around a Schwarzschild black 
hole with mass parameter M. The particle's worldline is denoted x£(r), where r is the proper time. The particle's 
trajectory is bounded within the range r m ; n < r < r max for certain r max and r m i n > 4M. Without loss of generality 
we shall take the orbit to be equatorial, namely 9 p — n/2. 

Assume now that the particle carries a scalar charge q. This charge couples to a massless, minimally-coupled scalar 
field <&(x M ), satisfying the field equation 

□$ = -47rp. (2) 
Here p is the scalar charge density, which takes the form of a <5-function along the particle's worldline: 

/oo 
6±[ x - Xp { T )][-g{x)]- l / 2 dT, (3) 
-oo 

where g = — r 4 sin 2 9 is the metric determinant, and hereafter the vectorial indices of x^{t) and x^(t) are suppressed 
for brevity. 

Since t is timelike (hence monotonic) we can use it instead of r to parametrize the orbit. In the r, t plane the orbit 
is then denoted by r = r p (t). Transforming the integration variable in Eq. © from r to t and substituting 9 p = tt/2, 
we find 

p = qir 2 ^)- 1 5[r - r p (t)] 5[<p - Vp (t)} 8[6 - n/2], (4) 

where u* = dt p /d,T. Note that it* only depends on r p (t): We have u* = E[l — 2M/r p (t)]~ 1 , where E is a constant of 
motion. 



B. Spherical-harmonics decomposition 

We now separate the field equation ^ by decomposing <f> in spherical harmonics, in the form 

oo I 

® = J2 12 <Pim(t,r)Y lm (d,tp)/r. (5) 



1=0 



Here Y\ m are the standard (complex- valued) normalized spherical harmonics given by Yi m — ci m Pi m (cos 9)e lmv , where 
Pi m are the associated Legendre polynomials and ci m are (real) normalization constants. The factor 1/r is introduced 
for later convenience. The charge density in Eq. ^ is decomposed in a similar manner: 

oo I 

P = ^2 12 pi m (t,r)Y hn (9,ip), (6) 

1=0 m=-l 

where 

(t, r)= f pY t * m dn = cb^rVrV^W^r - r p (t)]. (7) 



Here Q m = c; m Pz m (0), d£l = sin9d9dip is an element of solid angle, and an asterisk denotes complex conjugation. The 
field equation ([2]) is now separated, and for each Im the function 4>im(t,r) satisfies the partial differential equation 

d 2 4>i m d 2 (j>i m 

-7p Vl ( r > Vim = -47rr/ (r)p lm , (8) 
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where f(r) = 1 — 2M/r, r* = r + 2Mln[r/(2M) — 1] is the tortoise radial coordinate, and 

wm-/m(£ + ^). m 

The function <f>i m (t,r) is determined for each Zm by Eq. (|5J), supplemented with suitable boundary conditions at 
null infinity (no incoming waves) and at the event horizon (no outgoing waves) . Since the source term has a S- function 
form, the function (f>i m (t,r) must be continuous at r — r p (t), but it will generally fail to be differentiable there: Its 
r derivative (and also its t derivative, except at the two orbital turning points) will suffer a discontinuity along the 
orbit. On each side of the worldline r = r p (t), however, the function (j>i m {t, r) satisfies the homogeneous part of Eq. 
©, and, since the homogeneous field equation and the curve r — r v (t) are both analytic, we expect 4>i m {t,r) to be 
analytic (in both r and t) anywhere off the worldline. We shall assume this analyticity here, although we are not 
aware of a proof. The alternative option appears highly unlikely, because it would mean that the actual field produced 
by the point charge somehow develops irregularities in the vacuum region off the particle. 



C. Self force via mode sum 



Once the functions <pi m have been determined (e.g. numerically), the self force acting on the scalar charge may be 
constructed by the mode-sum method. This procedure is described in detail in Refs. [1, SOU. Here we outline it 
briefly, in order to provide some perspective on how the quantities 4>im are incorporated in the construction of the 
self force. 

Let denote the contribution of an individual I to $: 

l 

$i(t,r,e,<P)= Mt,r)Y lm (9,v)/r. (fO) 
m=— I 

The Z-component of the full-force field, Fin (x) , is then obtained by applying a certain linear differential operator 
to (pi. [J 7 ^ is the same differential operator which determines the force that a smooth, "non-self" field would 
exert on a test charge.] In the scalar-field case we have T n = qd u , and therefore [29[ 

l 

F hl = Fy&i = q Wi m (t,r)Yi m (B,<p)/r]. (If) 

rn= — I 

The quantities Fi^ are to be evaluated at the particle's location. To be more specific, let us denote by Xf the event 
(on the particle's worldline) at which the self force is to be evaluated. Then the quantities F^ are to be evaluated at 
x — > Xf (a somewhat rough statement which will be refined shortly). From Eq. (If I p it is obvious that Fi u depends 
linearly on the following functions of r and t (for each m): 4>i m , 4>i m ,n aud 4>i m ,t- We use the symbol <j)\ m (i — 0, 1, 2) 
as an abbreviated notation for these three key functions. 

As was discussed above, the quantities <fiim,r and (j>i m ,t are n °t truly defined on the curve r — r p (t), and in 
particular at x = Xf. Instead, each of these functions has two well-defined (but generally different) one-sided limits, 
corresponding to approaching the worldline point Xf from the range r > r p (t) or r < r p (t). We denote these two 
one-sided limits as i — > Xf+ and x — > respectively. Correspondingly, the quantities will each have two 

one-sided limits, Fjt and F,~ , defined by 



Ffc(x f )= lim Fi^x). (12) 



*Xf± 



The self force at x = Xf may now be derived from either set of quantities, Fjt or i<V , via the mode-sum formula Q 

oo 

°' f = E T (J + 1/2)4, - B»] , (13) 



1=0 



where and are certain parameters ( "regularization parameter" ) which Refs. [1,[T(3] determine analytically. Note 
that the two one-sided limits in Eq. (p~3j) yield the same value of F„ . 
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D. Frequency-domain analysis 



The partial differential equation ([8]), which determines the functions 4>i m (t,r), may be tackled in either the time 
domain or the frequency domain. In time-domain calculations, one directly integrates this equation numerically, 
using time-evolution on a two-dimensional grid. In the frequency- domain method, on the other hand, one first further 
separates this equation into Fourier frequency modes, using 



4>im{t,r) = I dto Ri muJ {r)e- wt (14) 

and 

-*irrf(r)fr m (t,r) = I duj Z lmuJ (r)e- wt . (15) 



Equation ((8|) then reduces to the ordinary differential equation 

d 2 Rl m L 



drl 



[Vj(r) - Lo']R lmuJ = Z lmuJ . (16) 



Since pi m {t,r) only has support on the curve r — r p (t), it follows that Zi muj (r) is only supported within the range 

^*min ^* '"max- 

From Eq. ([7]) it is evident that pi m only depends on t through r p (t) and (fp(t). For an eccentric geodesic r p (t) 
is periodic, and ip p (t) also has its inherent 2ir periodicity. It then follows (see App. [D]) that pi m is 2-periodic in t, 
namely it has a discrete spectrum of the form lo — nfl r + mCl v = Lu nm . Here fl r and Q, v are the two fundamental 
frequencies associated with the particle's radial and azimuthal motions, respectively. The integrals in Eqs. (fl"4"l) and 
(1151) thus reduce to summation over n. In particular, 

0j ro (t,r) =J2 R ^n m (r)e- iUnmt , (17) 

n 

where in principle the summation is over all integer values of n. 

The physically-acceptable solutions of Eq. (JTSJ) are those satisfying the appropriate boundary conditions at both 
edges r — > oo and r — > 2M, which correspond to pure outgoing waves at spatial infinity and pure incoming waves 
at the event horizon. The standard procedure for constructing the desired physical solution, for given Imu), begins 
with the construction of a basis of two independent homogeneous solutions, Rl muJ and R7 mw - These two homogeneous 
solutions satisfy the required boundary conditions at, respectively, r —> oo and r — » 2M (note that there is no 
non-trivial homogeneous solution which satisfies the required boundary conditions at both edges). One then utilizes 
the standard Wronskian-based formula for generating inhomogeneous solutions to second-order linear differential 
equations. Transforming the integration variable from r* to r using dr / dr* — f(r), and recalling the bounded support 
of Zi muJ {r), this formula takes the form 

R, (r) - R+ (r) f R i™^ r ') Z ^(r') , ( ) /~ W Rfmu,(r')Z lmu (r') , 

lJ ^ l) L Wf(r') dr+H im^ r )J r wf[rl) dr 

= Bfitir), (18) 

where 

W = R^ (dRLjdr*) - Rt mw (dR^Jdn) = const (19) 
is the Wronskian. In the regions r < r m j n and r > r max this formula reduces to the homogeneous solutions 

{^lmw^imu( r ) = ^ImuW' T — r min, 
(20) 
ChnujRlmLo( r ) — Rlmw( r )> r — r max, 

where the coefficients and Cfc are given by 

= w -l £~ RLJr)Z^(r) ^ (2l) 
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We conclude this section by explicitly writing the frequency-domain expressions for the three key functions <ft\ m , in 
the particle's neighborhood: 

<Plm(t,r) = ^2R^ Un Jr)e-^ \ (22) 



0lm,r(*,r) = Tr^- {r)e ~ %U)nmt > (23) 



4i m ,t(t, r) = -i^UnmR^Jr) . (24) 

n 

Since the particle resides in the range r m ; n < r < r max , the radial functions fi/Sf w (r) involved in these expressions are 
truly inhomogeneous [unlike the functions RxL^if) defined in Eq. (|20[) . which are homogeneous] . (30j 

III. THE HIGH-FREQUENCY PROBLEM 
A. Statement of the problem 

The functions (ft] m (t,r), whose one-sided values are required for the self- force calculation, are perfectly (one-sided) 
smooth, even at the coincidence limit r — > r p (t). Owing to this smoothness, one may naturally expect that there 
ought to be a way to calculate the required one-sided quantities in the frequency domain without referring to large cu 
values. Such a method indeed exists, as we explain in the next section. However, with a straightforward application 
of the standard frequency-domain method, based on Eqs. (|22p - (|2"4"|) . one finds that the Fourier series either fails to 
converge to the correct values (for 4>im,r and 4>im,t) or converges very slowly (for 4>i m ) as r — > r p (t). 

To demonstrate this convergence problem we consider first the Fourier sum (|23p for (fti m ,ri which is required for 
calculating F^ cU . [Essentially the same argument applies to Eq. (|24|) for <fti m ,t-] Suppose that we attempt to evaluate 
4>im,r at a point x = Xf on the worldline, with coordinates t = tf and r = Tf = r p (tf). The values of the radial 
functions dR^^/dr at r = r/ are just the Fourier components of the function 

<t>L, r (t) = 4> lm , r (r = r f ,t). (25) 

Since (j>i m ,r is discontinuous at the worldline, cft( m r (t) is discontinuous at t = tf [as well as at any other t value for 
which r p (t) = rj]. We therefore encounter here the Gibbs Phenomenon [2l|: If a function F(t) is discontinuous, its 
Fourier sum will fail to converge to the correct value at the discontinuity. (Away from the discontinuity the Fourier 
sum will converge, but rather slowly and only conditionally: The n-th order term will behave as 1/tj.) 

Consider next the convergence properties of the Fourier sum for tpi m in Eq. (f2"2"| . Defining <f>{ m (t) = 4>i m (r = rf,t), 
we observe that 4>{ m {t) is continuous, yet its derivative is discontinuous at t = tf. Standard Fourier theory [20] then 
has it that the Fourier sum will indeed converge (to the correct value) at r — > r p (t), but this convergence will be 
rather slow: The n-th term of the Fourier series is expected to behave as 1/n 2 . 

B. Numerical illustration: Scalar-field monopole 

It is instructive to illustrate the above problem with an explicit calculation. For this, we consider the example of 
the monopole mode I — m = 0. The spectrum in this case becomes simply w = nQ r (with integer n), and the Fourier 
sum (fTTf takes the form 

oo 

0(r,t)= £ CVje- 1 *'. (26) 

n— — oo 

We hereafter use the notation </>; =TO= o = (ft, Ri =m =o,uj=nn r = Rn, etc. to represent the various monopole quantities. 
For a given orbit, the inhomogeneous n-mode radial functions i?Jf h (r) can be computed numerically based on Eq. 
(If 8p . The relevant numerical procedure is rather standard, and we relegate its description to App.[C] In the following 
we present sample results and discuss their significance. 
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r/M r/M 

FIG. 1: Construction of the scalar-field monopole and its r derivative as a sum over inhomogeneous frequency modes in the standard 
approach. The numerical solutions shown here correspond to an eccentric geodesic orbit around a Schwarzschild black hole, with semi-latus 
rectum p = 10 and eccentricity e = 0.2 (see App. lAl. The "periastron" and "apastron" for this orbit are at r m ; n = (25/3)M = 8.333A/ 
and r max = 12. 5M, respectively. We used numerical integration to calculate the inhomogeneous radial functions H™"(r) through Eq. 
| |C4| |. and then obtained the partial sums <f>(t, r; n max ) and <t>,r(t, T\ n max ) as defined in Eqs. H27H and H28I I. Plotted here (solid lines) are 
the partial sums (per unit scalar charge) for n max = 2,4,8, 15, as functions of r at a fixed time t when the particle is at r = lOAf (this 
corresponds to radial phase \ = n/2; see App. [A}. The left panel displays the scalar field itself; the right panel shows its r derivative. For 
comparison, we also display (dashed line) the full scalar monopole solution, which we obtained directly using numerical evolution in the 
time domain (for our purpose, this latter solution can be taken as an accurate benchmark). It is evident that the n-mode sum converges 
quickly to the correct value in the regions r < r m j n and r > r max , but the convergence deteriorates inside the domain r m ; n < r < r max , 
where "Gibbs waves" set in. The partial sum over frequency modes is smooth at the particle, and hence, strictly speaking, cannot recover 
the true jump discontinuity in the field derivative there. Finding a way around this technical problem is the main goal of this work. 



Figures [T]and [5] display numerical solutions for the sample orbital parameters r max = 12. 5M and r m i n = (25/3)M = 
8.333M. (This corresponds to "semi-latus rectum" p = 10M and "eccentricity" e = 0.2, both quantities defined in 
App. [All In Fig. [T] we plot the (real- valued) partial sums 

n max 

<t>{r,t;n m ^) = ]T iC h (r)e-™ n «* (27) 



and 



"max , 

Mr,*;n max )= £ -iCVK" 1 ""' (28) 



n— - 



as functions of r at the fixed time t when the particle is located at r — 10M, for a sample of n max values. For 
comparison, we also plot the full monopole solution (and its r derivative), which we obtain using a time-domain 
numerical evolution code similar to that developed by Haas in Ref. [16|. Evidently (and as expected), the convergence 
of the rt-mode sum for both <j> and cf> }r seems very fast at r < r m j n and r > r max , but deteriorates significantly in the 
domain r m j n < r < r max , where "Gibbs waves" dominate the behavior. Figure [2] illustrates the convergence properties 
of the partial sums </>(n max ) and ir («max) at the very location of the particle (on the same time slice as in Fig. [T]). 
The data on the left panel suggest that the partial sum for the field <fi converges at the particle as ~ l/n max — in 
accordance with theoretical expectation [20]. The results shown in the right panel of Fig. [2] demonstrate that the 
partial sum for the derivative <p tr fails to converge to the correct (one-sided) values. They also (loosely) suggest that 
this partial sum in fact con verg es to the two-side average value of 4>,r at the particle. This, indeed, would again accord 
with theoretical prediction [20j. 



C. Practical implications of the high-frequency problem 

The above numerical example serves to illustrate the following: From a practical point of view, it would seem very 
difficult to extract the correct values of the key functions <p\ at the particle's location, based plainly on a naive 
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left-hand values 




max max 



FIG. 2: Convergence of the Fourier re-mode sum at the location of the particle (illustrated here based on the numerical data of Fig. [TJ. 
Left panel: The deviation (per unit scalar charge) of the partial sum </>(n max ) from the full field <j> at the location of the particle, as a 
function of n max . The n-mode partial sum appears to converge to the correct value approximately as l/n max (note the log-log scale), 
as expected on theoretical grounds. Right panel: The deviation (per unit scalar charge) of ^ r (n ml ,) from the full-field derivative r , 
again evaluated at the particle's location. Since the true <f> >r [unlike ^ r (?in«)] has two different one-sided values at r = r p , so does the 
deviation. The two solid lines represent these two one-sided values of the deviation. (These two values actually have opposite signs for 
each n max , which is obscured here since only the absolute value of the deviation is shown.) From the essentially- horizontal shape of the 
two solid lines it is evident that the n-mode sum for A f does not converge to the correct one-sided values. As an aside, we also plot here 
(dashed line) the difference between </>,r-(n max ) and the two-side average value of the full derivative <j> r (per unit scalar charge, at the 
particle's location). The graph loosely suggests (referring to its upper envelop and ignoring the seemingly oscillatory deep structure) that 
<^>,r(fi max ) converges, albeit very slowly, to the average value of <j> lT at the discontinuity. That, indeed, would again be consistent with 
theoretical prediction. 

summation over frequency modes as in Eqs. The partial Fourier sum for the field 4>i m itself would converge 

very slowly (as 1/ro) and its evaluation would hence be computationally expensive. Worse, the partial sums for the 
derivatives 4>i m ,r and (f>i m ,t would simply fail to yield the desired one-sided values at the particle, even if one could 
sum over infinitely many modes. 

Having stated the above, we should also point out that Eqs. and (f!M]) should not be deemed entirely useless 
for the purpose of calculating 4>i m , r and (t>i m> t at the particle. In principle, one could pick a point close to r = r p (t), 
yet not quite at r p (t), and calculate 4>im,r (say) there. The Fourier sum will converge at this point, although rather 
slowly. One could then pick a series of r-values which approach r p (t) (say, from the '+' side), calculate <f>i m ,r at each 
of these points, and then evaluate the desired sided-limit of 4>im.r through extrapolation. This procedure, however, is 
cumbersome and is hardly likely to be computationally tractable. 

An alternative implementation strategy could make use of the fact that the sums in Eqs. (J23]) and actually 
converge to the two-side averages of (f>i m ,r and <pi mi t at the particle (recall the discussion relating to Fig. [5]). These 
average values (along with <f>i m itself) could be used to construct the average force modes = (F^ + Fj~)/2, 
which could then be directly implemented in a "two-side averaged" version of the mode sum formula: i* 1 ® = 
Yli^o — Although this method is likely to be by far more efficient than the extrapolation method 

mentioned earlier, it would still present a computational challenge, as this method, too, involves the evaluation of 
slowly-converging Fourier sums. 

IV. METHOD OF EXTENDED HOMOGENEOUS SOLUTIONS 

A. Formulation of method 

In this section we describe our alternative method for frequency-domain construction of the key quantities (j)\ m 
required for calculation of the self force. This method, to which we refer as the method of extended homogeneous 
solutions, completely avoids the high-frequency problem described above and ensures exponentially-fast convergence 
of the Fourier series. 
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We begin by extending the definition of the homogeneous functions Ri m(jj { r ) m Eq. (120]) to the entire domain 
r > 2M: 

£L>) = CjL*L>). r>2M. (29) 

We then define the two time- domain extended homogeneous solutions 0^ and <f)J to be the outcome of replacing 
in Eq. |22|) by the homogeneous solutions ii^, or Rf muj , respectively: 

~^r)^Y. R t^Jr)e—\ (30) 

n 

We emphasize that each of the fields <j>i~ and <j>J is defined in the entire domain 2M < r < oo. 

The convergence properties of the sum in Eq. (|30f> are dictated by the large- |n| asymptotic behavior of the coefficients 
-R; mw (r). This high-frequency asymptotic behavior can be examined using a WKB-type analysis, which we carry 
out in App.|Dj This analysis shows that, at least within the leading-order WKB approximation, the terms on the right- 
hand side of Eq. (|30[) decay (at least) exponentially in \n\. This exponential decay is uniform in t and r, throughout 
r > 2M. Also, one naturally expects that the contribution from higher-order terms in the large-w WKB expansion 
will converge even faster than the leading-order contribution, and hence will not affect this uniform exponential decay. 

The exponential convergence of the sum (|30[) is extremely convenient for numerical applications, as illustrated in 
the next subsection. But it also has important mathematical consequences. Since the homogeneous radial functions 
flit, (r) are analytic, the uniform exponential decay of the individual terms in the above sum implies that the overall 
sum — namely the extended homogeneous solution — is an analytic function of r and t throughout r > 2M. 

We now argue that on each side of the curve r = r p (t) the actual time-domain function 4>i m (t, r) coincides with one 
of these extended homogeneous solutions, namely 

f r>r p (t), 
<hm(t,r) = l (31) 
Ui^(*,r), r<r p (t). 

For concreteness, let us present our argument explicitly referring to the first of these equalities: (i) In the domain 
t > i"max this equality obviously holds because itjjj£, and Rf mu coincide in that domain, (ii) As was already mentioned 
in Sec. HH we assume that the function <fii m {t,r) is analytic throughout the range r > r p (t) [as well as in the other 
range, 2M < r < r p (t)], because the alternative option appears unreasonable, (iii) As was just discussed above, the 
high-frequency analysis in App. [D] strongly suggests that the extended homogeneous functions </>/L (£)?*) are analytic 
throughout r > 2M. (iv) Since both functions <j>i m and <^ m are analytic throughout the domain r > r p (t) [from (ii) 
and (iii)], and coincide at r > r max [from (i)], they must coincide throughout r > r p (t). (v) By continuity of both 
functions tj>i m and <^L, they coincide at r — r p (t) as well. Obviously, the same line of argument applies to the second 
of the equalities (|3"Tj) as well. 

In the rest of this subsection we describe the utility of the extended homogeneous fields defined above in calculations 
of the self-force via the mode-sum method. Recall from Eqs. (fTTji - fTB)) that this method requires as input (either of) 
the one-sided limits of 4>\ m at the particle, which we now denote 

<t>t(x f )= lim 4 m . (32) 

x — yx f ± 

In terms of <£J^, the various components of the quantities Ff^(xf) [as defined through Eqs. (fTTj) and (fT2|) and used in 
the mode-sum formula fTS")) ] are expressed directly as 

i 

= £ E {^m,V 4L,r ^ 4>L/rf> im 4m} Yi m ^/2,<p f ) (33) 

' m=-l 

(along with Fj^ = 0). 

For concreteness, let us focus first on one of the quantities say (f>f m - By definition, the limit x — > Xf+ in Eq. 
(|32l) only samples the range r > r p (t). Using Eqs. (|3T|) and (|30|) we may re-express 4>f m as 

4>Uxf) = Hm E^U>)e- JlW . (34) 

n 
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Since the sum over n here converges uniformly, we may interchange the limit and summation. However, the functions 
Rf muJ and e~ luJnmt are analytic, so we can now omit the limit x — > xj+ and instead simply evaluate these functions 
at x = xj. The final outcome from these manipulations is stated in (the '+' case of) Eq. ([55)) below. 

The above treatment is equally applicable to <t>~[ m , and we obtain a similar formula for constructing ^^{xf) out of 
the extended homogeneous modes (r) [the ' — ' case of Eq. ([55]) ]. Moreover, the same treatment also applies to 

film r an d ^im v The six key quantities (j)]^ can all be constructed from the extended homogeneous radial functions 
(and their derivatives) in the form 

O/) = E^J^e---'/, (35) 

n 

tfmA x f) = E ^ r R L^Jr f )e- l — tf , (36) 

n 

O/) = - i T / u ^L, nm (Tf)^ m * f - ( 37 ) 

n 

Equations (p?5|) and ([33)1 - (|3T|) . combined with Eqs. ([3"3")l and (fP3")l , constitute our new method of calculating the self 
force in the frequency domain. The high-frequency problem is entirely circumvented in this method, as the Fourier 
sum converges exponentially- fast for all functions cftVz. 



B. Numerical illustration: Scalar-field monopole revisited 



Let us revisit the calculation of the scalar-field monopole — this time using the method of extended homogeneous 
solutions. The homogeneous basis solutions K„(r) are constructed numerically in just the same manner as in the 
standard approach (see App.jBj. Then, however, instead of calculating the actual inhomogeneous modes R™ h (r) as m 
Sec. IIII1 we construct the extended homogeneous solutions R„[r) as they are defined in Eq. (p?9")) . with the coefficients 
Cz= m =o uj=n,n = calculated through Eq. (jC7[) of App. [C] The time-domain extended fields and their r derivatives 
are then approximated by the (real-valued) partial sums 

Tl mas 

± (r,t;Tw)= E Rt(r)e- inQnt , (38) 



"max j 

<^(r,t;n roax ) = E ^M^" 10 " 4 , (39) 

"— — "max 

with sufficiently large n max . 

We point out the following matters relating to the implementation of Eqs. ([3"B")) and ([3"9"| : (i) In the new approach, 
the computation of <f> (or (f> jT ) for all r and t involves the (numerical) evaluation of only two integrals — the ones in Eq. 
(|C7[) of App. [Cj In contrast, the standard approach requires the evaluation of two integrals — the ones in Eq. (|C4[) — 
separately for each value of r between r mm and r max . (ii) The full scalar monopole is continuous at r = r p (t); however, 
the contributions to the extended functions 4» + and (j>~ from each individual n mode do not match continuously along 
this curve. Consequently, for any finite rt ma xj the partial sums for these extended functions do not match at r = r p (t). 
The amplitude of this mismatch is expected to decrease rapidly with growing rtmax, as the results below indeed 
demonstrate. 

Figures [3H5] display numerical solutions obtained based on Eqs. ([3"8")) and ([3"9"| . Our goal here is to assess the 
performance of the new method against the standard method, and to this end we have chosen for our numerical 
experiment the same orbital parameters as in Figs. [1] and [D of Sec. Mil For clarity, we only show the '+' and ' — ' fields 
in their respective relevant domains, i.e., r > r p (i) for the former and r < r p (t) for the latter. 

Our numerical illustration serves to demonstrate the following: (i) The sum over '+' and ' — ' extended n-modes 
converges quickly to the correct, full solution everywhere in the respective domains r > r p (t) and r < r p (t). (ii) In 
particular, the mismatch between the values of the '+' and ' — ' partial sums at the particle's location quickly converges 
to zero with growing rt max . (iii) The convergence of the extended n-mode sum is exponential everywhere — even in the 
region r m j n < r < 7" max ; in particular, it is exponential at the very location of the particle. This applies to both the 
field and its derivatives, (iv) The new scheme completely circumvents the Gibbs effect which disrupts the convergence 
of the inhomogeneous n-modes in the standard approach. 
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FIG. 3: (To be compared with Fig. [T]) Construction of the scalar-field monopole (left panel) and its r derivative (right panel) using 
"extended homogeneous solutions". The orbital parameters are chosen here as in Fig. [T] i.e., p = 10 and e = 0.2, and we again present 
the various fields as functions of r at a fixed time t when the particle is at r = lOAf. Vertical dotted lines mark the particle's periastron 
and apastron radii, r m ; n = (25/3)M and r max = 12. 5M, respectively. The solid line represents the full scalar monopole solution, obtained 
using numerical evolution in the time domain. The broken lines represent partial sums over extended homogeneous n-modes, calculated 
numerically based on Eqs. Q38|l and l|39|l . For clarity, in both panels we plot the '+' and ' — ' partial sums only in their relevant domains, 
r > 10M and r < 10M, respectively. We show here the partial sums for n max = 0, 1, 2 only — the partial sums <^> (n ma x = 3) are already 
indistinguishable from the full solution at the scale of this plot (but see Fig. [4] below) . The individual n-modes of the extended fields 0+ 
and <j>~ do not match continuously at the location of the particle, but their sum seems to converge quickly, everywhere, to the true solution 
(which is continuous). Similar fast convergence is manifest also for the derivative >r . "Gibbs waves", which disrupt the convergence of 
the actual inhomogeneous n-modes, are altogether avoided within the new scheme. 




FIG. 4: Convergence of the extended n-mode sum for the fields (j)^ (left panel) and their derivatives j7 - (right panel). For the same 
case shown in Fig. [3] we plot here the fractional differences between the partial sums and the full field (or the full-field derivative), for 
n max = 1-5. The various graphs are labeled by their corresponding n max values. The middle vertical dotted line marks the particle's 
momentary radius at r = 10M. Once again, we display the '+' and ' — ' values only in their respective relevant domains r > 10M or 
r < 10M. Note the exponential scale of the y-axis. (The seemingly odd behavior of the data for n max = 2 and n max = 5 in the right panel 
is simply due to a change-of-sign which the corresponding fractional differences happen to experience around r = 1AM and r = 10M, 
respectively. The tiny wiggly feature, barely visible near r = 10M for n max = 5, is due to the numerical error in the time-domain 
data, which for <t> >r is estimated at ~ 10 — 6 in fractional terms.) The exponentially-fast convergence of the extended n-mode sum even at 
frnax — and particularly at the very location of the particle — is evident from these plots. 
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FIG. 5: Convergence of the extended n-mode sum at the particle's location. The left and right panels display the values at the particle's 
location (r = 10M) of the fractional differences shown in the left and right panels of Fig. U correspondingly. The left panel demonstrates 
the exponential convergence of the extended n-mode sum to the correct value of the field at the particle. In particular, the mismatch 
between the partial sums for 0+ and <f>~ at the particle appears to vanish exponentially with increasing n max . As we suggest in Sec. 
FVl one can in fact make a good use of this mismatch in numerical calculations, by recording its residual value and using it to assess the 
truncation error of the partial n-mode sum. Comparing with Figs. [T] and [2] (left panels), it is striking that, in the example considered here, 
summing up to only n = 2 with the new method achieves a better accuracy in the local field than summing over as many as 16 modes 
using the standard approach. The right panel demonstrates the fast convergence of the derivatives and cf>~ r in the new approach. An 
exponential convergence, expected form theory, is loosely suggested from the data shown. 

V. DISCUSSION 

The basic relation underpinning our new method is expressed in Eq. (|3ip . which describes the construction of the 
Zm-multipole of the scalar held from extended homogeneous solutions. Eqsuation (|29|) and (|35| - ((37|) . in conjunction 
with Eqs. (|33[) and (|13[) . constitute our new prescription for constructing the self force in the frequency domain. The 
following list highlights the advantages of this new formulation. 

• In the standard frequency-domain scheme, the Fourier sum over w-modes suffers from Gibbs-type irregularities 
near the particle's location. In particular, the Fourier sum fails to correctly recover the derivatives of the field's 
multipolcs at the particle (which are essential input in self- force calculations). In the new method one constructs 
the local field's multipoles as Fourier sums of globally homogeneous aj-modes. These sums converge uniformly, 
circumventing the above complication. This is the essential and most crucial merit of the new method. [Note 
also that the uniform convergence (which also applies to the derivatives of the extended homogeneous mode 
functions) allows one to obtain the derivatives of the field's multipoles using a term-by-term differentiation of 
the individual Fourier components. This again leads to Eqs. (|3l))) and (|37|) above.] 

• The sum over the homogeneous lo modes converges exponentially everywhere, and even at the very location of 
the particle. This is extremely convenient from a practical point of view. It should be noted that, within our 
new scheme, it becomes "as easy" to obtain the local /-mode field near the particle as it is to obtain the same 
field in the far-zone — in sharp contrast with the situation in both the standard frequency-domain method and 
the time-domain method. 

• In the standard scheme, each of the inhomogcneous w-modes is computed via Eq. (|18[) . In practice, this involves 
the (numerical) evaluation of two integrals for each value of r between r m i n and r max . In the new scheme, the 
extended homogeneous w-modes are obtained via Eq. (|29[) . which requires the same two integrals (|21[) for all 
values of r. Thus, remarkably, the new scheme does not only perform better mathematically — it is also much 
simpler to implement. 

• Finally, within the new scheme one is offered a convenient handle with which to monitor and control the large-w 
truncation error (i.e., the error caused by omission of the terms |n| > n max ): The residual amount by which the 
partial sum (f° r an y °f the quantities <fii m ) fails to be consistent with the appropriate jump condition 
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at r = r p (t) is a faithful measure of the truncation error. One can therefore conveniently keep the latter below 
a set level by setting a threshold on the amount of residual inconsistency. 

What is the application scope of the new method? In this paper we introduced the method as applied specifically 
for scalar-field perturbations on the Schwarzschild background. However, the basic idea is equally applicable to a wide 
range of other problems. It is clear from our discussion that the problematics of reconstructing the local multipole 
field as a sum over frequency modes has little bearing on the precise form of the underlying field equations. Rather, 
it is a feature of the spherical-harmonic decomposition: When we consider an individual Im mode, we effectively 
convert the physical problem of a point particle moving in an eccentric orbit into a problem of a source "shell" (which 
expands and contracts over time); the perturbation field is not smooth across the shell, which gives rise to Gibbs-type 
complications when we attempt to express it as a sum over frequency modes. The same problem would occur in 
essentially any perturbation treatment in Schwarzschild which incorporates a spherical harmonic decomposition, such 
as the Teukolsky formalism (for EM and gravitational perturbations), the standard Regge-Wheeler/Zerilli/Moncrief 
formalism of gravitational perturbations, or the more direct Lorenz-gauge formulation [231 ]. The method we propose 
here as a cure for the problem is directly applicable for any of these treatments. A forthcoming paper [191 ] will 
report on the calculation of the monopole and dipole modes of the Lorenz-gauge metric perturbation (for eccentric 
orbits in Schwarzschild), facilitated by the method of extended homogeneous solutions. This calculation is now being 
incorporated in a code which computes the total gravitational self- force for generic orbits in Schwarzschild j24| . 

To what extent is our new method relevant for Kerr perturbations? Here the situation is more subtle. The original 
mode-sum scheme for the self force [lfj| (which sets the main context for the current work) incorporates a decomposition 
in spherical harmonics even in the Kerr spacetime. This is, of course, technically possible (see the footnoted remark 
at the Introduction), although the resulting field equations then couple between different i-modes. Regardless of the 
latter fact, each of the Zm-modes in this decomposition would again be sourced by an expanding/contracting thin 
shell, the non-smoothness of the perturbation across this shell would give rise to the Gibbs phenomenon, and our 
method would provide an efficient cure. 

The situation is different if one tackles the Kerr problem by means of the more natural decomposition in spheroidal 
harmonics, which decouples the field equation in the frequency domain. Since the spheroidal-harmonic functions 
depend on the frequency, one no longer has a strict notion of a time-domain L lm mode' in this case. One may 
(somewhat artificially) define the "spheroidal-harmonic /'m-mode" ^i' m (t, r, 8, if) by summing over all u> for given 
spheroidal-harmonic numbers l',m. However, in this case the effective geometric picture of a thin source shell would 
no longer apply: In the procedure of Fourier decomposition of the original point source (to obtain the source's I'mu 
modes) followed by a Fourier summation over to, the extra dependence of the spheroidal harmonics on lo will cause 
the reconstructed source's I'm mode function to deviate from a 5- function in t (for given r, 6, ip). Correspondingly, at 
a given t the source's I'm mode function will most likely represent a "smeared" shell. 

Based on the above discussion one might conclude that the high-frequency problem would not occur in the first place 
if spheroidal harmonics were used (as in this case there would be no <5-type shell). We believe, however, that this may 
represent a false logic. The viability of the mode-sum approach relies crucially on the fact that the individual mode 
contributions Fjj t in Eq. (fT3"|) are well-defined quantities. This fact is a direct consequence of the perfect one-sided 
smoothness of the mode functions <pi m {t,f) even at the limit r — > r p (t). This smoothness, in turn, stems from the 
fact that for each spherical- harmonic /, m the source term is confined to a S- function over a shell — and this 5- function 
is distributed over the shell in a perfectly smooth manner. Note also that the functions 4>i m {t,r) are homogeneous 
time-domain solutions on both sides of this shell — even at the immediate particle's neighborhood. In spheroidal- 
harmonic decomposition this situation is changed, and the spheroidal-harmonic mode function $;< m (t, r, 9, ip) defined 
above will no longer be a homogeneous solution in the very neighborhood of the particle. It is conceivable that in 
the immediate particle's neighborhood the smeared source shell will be dominated by large-w modes. [The spherical- 
harmonics decomposition is protected against this potential problem thanks to the combination of (i) the perfect 
off-shell homogeneity, and (ii) the independence of the I, m harmonic — and hence of the ^-function distribution over 
the spherical shell — on to.] Thus, in a spheroidal-harmonic decomposition, the high-frequency problem may take a 
much more severe form: It may endanger the very existence of (namely the spheroidal-harmonics analogs of Ff^) 
as regular quantities, which would in turn render the (spheroidal-harmonics analog of the) mode-sum formula (|13|) 
meaningless. 

The morphology of the spheroidal-harmonics smeared source shell still needs be investigated, and especially its 
structure near the point charge. The mode contributions F^f may turn out to be well-defined after all, but this is 
far from obvious. In any case, a spheroidal-harmonics-based variant of the mode-sum method for Kerr has not been 
developed yet to the best of our knowledge. The existing Kerr mode-sum variant [10j incorporates the spherical- 
harmonics decomposition, and as such it exhibits the same high-frequency problem as in the Schwarzschild case. The 
method of extended homogeneous solutions then elegantly resolves this problem in the Kerr case as well. 
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APPENDIX A: ECCENTRIC GEODESICS IN SCHWARZSCHILD 

This appendix reviews the standard description of eccentric geodesies in Schwarzschild spacetime, and provides the 
necessary formulas for all orbital quantities (frequency, four-velocity, etc.) needed for the numerical computations in 
Sees. ED and [TV] 

As in the main text, we consider a pointlike test particle in a bound equatorial geodesic orbit around a Schwarzschild 
black hole with mass parameter M. The radial location of the particle is bounded in the range r m j n < r < r max , for 
some r max > r m i n > 4M. Such geodesies constitute a two-parameter family. Each geodesic is uniquely characterized, 
for example, by the two "turning point" values r m i n and r mayi . An alternative parameterization (originally due to 
Darwin [2a |) employs the "semi-latus rectum", p, and "eccentricity", e, both analogous to their counterparts from 
Keplerian celestial mechanics. The parameter pairs (p, e) and (r m i n ,r max ) are related through 



or, inverting, 



2^max^*min ^max ^"min / a -i \ 

P = ; 1 e = ; > ( A1 ) 



(A2) 



■ ihcia -1 5 ' 111111 -i , 

1 — e 1 + e 

With the parameterization (p, e), the orbital radius is given by 

r P (x) = ~\~r~~~ > (A3) 

1 + e cos x 

where \ is a monotonically-increasing parameter ("radial phase") along the worldline. This parameter is related to 
the Schwarzschild time t p along the worldline through 



dx _ (p - 2M - 2Mecosx)(l + ecosx) 2 (p/M - 6 - 2ecosx 
~dT v = p~ 2 ^ (p/M - 2) 2 - 4e 2 



1/2 

(A4) 



with the constant of integration fixed such that x = at some "periastron" passage (r p = r m ; n ). x is related to the 
proper time along the eccentric geodesic through 

dx _ (1 + ecosx) 2 ( p/M - 6 - 2ecosxV /2 ^ A5 ^ 



dr M(p/M) 3 / 2 \ p/M 
The radius r p is manifestly periodic in \i with t-period 



T r = / {dx/dt^dx (A6) 



o 



' max 



and radial frequency Q r = 2tt /T r . 

For our numerical implementation of Eqs. (|C4|) and (|C7[) (in App. [Cj , we start by choosing r n 
use Eq. (|A1|) to determine p and e [or, alternatively, we pick p and e, then use Eq. (|A2[) to determine r m i n and r n 
We then solve for t p (x) for < x < t by integrating the inverse of Eq. (|A4p numerically [taking t p (0) = 0]. We hence 
obtain the radial period T r — 2t p (x — 7r) and the radial frequency fl r . Using Eq. (|A3|) to express x m terms of r 
along the orbit (again for < x < k), we then obtain t p (x{r)) = i P (?'), and, inverting in the range < t p < T r /2, 
also r p (t). Finally, u l is obtained as a function of t by writing u* = (dt p / dx){dx/ dr), substituting from Eqs. (|A4|) and 
(|A5|) . and then using Eq. (|A3|1 to express x m terms of r p (t). This procedure yields all necessary orbital parameters 
and functions for our numerical examples. 

Since the numerical illustrations of this work focus on the monopole mode, which is axially-symmetric, they do 
not require an explicit computation of the azimuthal frequency For completeness, though, we mention that this 
frequency is defined, in an orbit-average manner, by 

1 " T - 



n v = — (dtp p /dt)dt. (A7) 

J r JO 
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The "local" frequency, which depends on r p , is given in terms of the r-phase x as 

difip ip/M — 2 — 2ecosx)(l + ecos%) 2 
~~df = M(p/M)V 2 [(p/M - 2) 2 - 4e 2 ] 1 /2 ' 

The frequency Q v can be computed by changing the integration variable in Eq. (|A7[) to x> using Eq. (IA4I) . 



(A8) 



APPENDIX B: NUMERICAL CONSTRUCTION OF THE HOMOGENEOUS SOLUTIONS 7?±(r) 

We describe here the numerical construction of the frequency-domain homogeneous basis (R^ , R„ ) m the example 
of the scalar-field monopole. For I = 0, the homogeneous part of the n-mode radial ODE (fT6|) takes the form 



^ + [n 2 nl - 2M/(r)r- 3 ] R$ = 0. (Bl) 

For each n, i?+ and R~ are independent solutions of this equation, satisfying physical boundary conditions at r — > oo 
and r — > 2M, respectively. These homogeneous solutions are needed both in constructing the 'true' inhomogeneous 
n-mode solution R™ h [via Eq. (|C4[) in App. [C] and in constructing the extended homogeneous solutions i?„ within 
our new method [via Eqs. (55]) and (jCT]) ]. 

Consider first the static mode, n = 0, which can be solved for analytically: We have, simply, Rq — r and Rq = r In /, 
which constitutes a unique basis (up to multiplicative constants) with the property that Rq is regular at the event 
horizon and Rq is regular at r — ► oo. [To see this, recall that the static mode of the actual, full scalar field is 
<f> oc Rq /r, with corresponding "internal" solution $~ oc const(^ 0) and "external" solution $ + oc In/. The former 
is regular at the horizon but fails to vanish at infinity, while the latter falls off as 1/r at infinity but diverges at the 
horizon.] 

For each mode n ^ we solve Eq. (|B1[) numerically with suitable boundary conditions, as we now describe. Our 
numerical domain is a one-dimensional array representing physical radii in the range r m < r < r out . The boundaries 
are taken to lie in the asymptotic vacuum domains: n n /(2M) -1 « 1 (the event horizon) and r ou t 3> M (spatial 
infinity). In practice, it proved sufficient for our purpose to set r; n = 2.001A/ and r out = 1000M. 

Consider first R~ . As an inner boundary condition for this function we use the ansatz 

R~(r) = e- inn ^ Y,<k(r-m k (at r = r ln ), (B2) 

fc=0 

where are coefficients to be determined below, and A:~ ax is taken large enough to guarantee that truncation error 
is kept below a prescribed threshold (in practice, fc~ ax = 10 proved sufficient for the values of e,p, n considered in this 
work). The oscillatory factor in Eq. (|B2[) is chosen such that the contribution from each n-mode to the full monopole 
field attains the asymptotic form oc exp[— inCl r (t + r*)] as r — > 2M [recall Eq. I|26p]. This represents purely ingoing 
radiation, which is the correct physical condition at the event horizon. To determine the coefficients a~ k we substitute 
Eq. (|B2[) in the field equation (|Bip and solve the resulting hierarchy of algebraic equations at each order in r — 2M. 
This yields the following recursion formula for the coefficients a~ k (with given n; we omit here the index n for brevity): 

x {[M{2k-l)(k-2)-l2iM 2 ntt r (k-l)}a k -_ l 



k>0 2M 2 k(k - UMnQ r ) 

+ [[k - 3)/2 - 6iMnO r ] (jfe - 2)a' k _ 2 - inn r (k - 3)a A T_ 3 } , (B3) 

with a k<0 — 0. All coefficients a k>0 are constructed recursively given a^, and are all proportional to Oq~ . 

To solve Eq. (|B1[) for R~(r), we simply set = 1, impose the values of R~ and dR~/dr at r- m using Eq. (|B2[) . and 
integrate numerically forward from r = r- m to r = r max (the value of R~ at r > r max is not needed in our analysis). 

Now consider R+. As an outer boundary condition for this function we take 

R+(r) = £ a + nk r~ k (at r = r out ), (B4) 

fe=0 

where a^ k are determined below, and &+ ax is chosen, once again, such that truncation error is kept below a set 
threshold (here, too, fc^^ = 10 was sufficient in our analysis). With this condition, the contribution from each n- 
mode to the full monopole field has the asymptotic form cx exp[— inQ r (t — r*)] as r — > oo, representing purely outgoing 



16 



radiation — the correct physical condition at infinity. The expansion coefficients in Eq. (IB4|) are obtained recursively 
using 

= 2^fc 1)a ^ + 2M{k 1?a+k -^ ' (B5) 

with a^ <0 — 0. (Recall n = is dealt with analytically, so one need not worry about the ill-definiteness of the 
recursion relation in this case.) All coefficients a^ >0 are constructed recursively given Oq , and are all proportional to 
a+. 

To solve Eq. (|B1[) for we set cl~q = 1, impose the values of i?+ and dR^/dr at r out using Eq. (|B4[) . and integrate 
numerically backward from r = r out to r = r m - ln (the value of at r < r min is not needed in our analysis). 



APPENDIX C: CONSTRUCTION OF THE INHOMOGENEOUS SOLUTIONS i?' n 



This Appendix details the computation of the inhomogeneous radial functions i?Jf h (r) for our numerical illustration 
in Sec. |HI1] 

For I = 0, the functions R™ h (r) satisfy the ODE 



r Yl pinh 

d -^- + [n 2 fi? - 2M/(r)r- 3 ] = Z n , 



(CI) 



where the source term is obtained by taking the inverse-Fourier transform in (the discrete, I = version of) Eq. {15 



Z n (r) = T- 1 f ' [-47rrf(r)p l=m=0 (t,r)}, 
Jo 



T.n T t dt 



g(4 7 r) 1 / 2 T- 1 / f (r^ru 1 )- 1 5[r - r p (t)}e inn ^dt. 



(C2) 

Here T r is the radial period [see Eq. (|A6j) ]. and in the second equality we have substituted for p from Eq. (O, 
setting m = and c/ =m= o = (47r) -1 / 2 . Note in Eq. ()C2[) that Z n = for r < r m ; n or r > r max , and that for any 
^min < t < fmax the integrand is supported only at two points within the integration domain — the two times t for 
which r p (i) = r. Changing the integration variable from t to r p one thus obtains 



Zn(r) = - 2q ^} [/^ cos[nn n t p (r)] x 9(r - r min ) x 6(r n 



r), 



(C3) 



where O is the standard unit step function, u r = dr p /dr is the r component of the particle's four-velocity, and t p (r) 
is the result of inverting r — r p (t) in the domain < t < T r /2, assuming r p (0) = r m i n [note r p (t) is single- valued in 
this restricted domain]. 

In the standard approach, the physical solution to Eq. (|C1|) is obtained through the formula (fTS]) , replacing Ri mul — » 
i?Jf h , -RrL — > -R^, and Z; mu — > Z n . Here i?+ and R~ are two independent solutions to the homogeneous part of 
Eq. (|C1|) . satisfying physical boundary conditions at r — > oo and r = 2M, respectively. It is convenient to change 
the integration variable in Eq. (|18p from r to t p (r) [taking i p (r m i n ) = 0], which avoids the singularity in Z n (r) at 
r = r max ,r min . Substituting for Z n from Eq. (| C 3 [) . and using |u r | = w r = u'(dr p /di) for < t < T r /2, Eq. (Tl8|) 
becomes 

r*»« H-(r p (t)) 



^CV) = -2 g (47r) 1/2 T,T 1 W r - 1 x 



+ i?-(r) 



o r p (t)u*(r p (t)) 
T - /2 fl+(r P (*)) 

tp(r) r p (<K( r p(*)) 



cos(ntt r t)dt 
cos(nfl r t)dt 



(C4) 



where we have introduced 



0, r < r min , 

*p( r ) = 4 *p( r )' r min < r < r n 



(C5) 
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The solution R™ h in Eq. (|C4|1 is manifestly an analytic function of r everywhere, except at r = r m i n ,r max . Note, 
recalling Eq. (|20[) . that outside the domain r m j n < r < r max Eq. (|C4|) reduces to the homogeneous solutions 

f Cn^imcW) r < r min , 
K lh (r) = < + + (C6) 

where the coefficients are given by 

Ct = -M^T-^ I"' 2 rp ^ff ( l )} cos (nn,*)*. (C7) 

Equation (|C4|) can be implemented numerically to obtain R^ h (r), in the following manner: We specify the physical 
orbit by picking the values of r m j n and r max , and then use the relations given in App. \X\ to obtain (numerically) 
the values of Q r and T r and the functions r p (t), t p (r) and u'(r p (t)) for the specified orbit. We next construct 
the homogeneous basis (r) by numerically integrating the homogeneous part of Eq. (jCip with suitable boundary 
conditions. This procedure is described in App. [Bl Once the solutions -R„(r) are at hand, the (constant) value of the 
Wronskian is obtained using Eq. (jTSJ) . Finally, for each given n, we calculate the integrals in Eq. (|C4[) numerically, 
and construct the solution R^ h (r). 



APPENDIX D: HIGH-FREQUENCY ANALYSIS 



In this appendix we analyze the convergence and analyticity of the Fourier sum involved in the construction of the 

(i) 

key functions (pi' in the "extended homogeneous solutions" approach — i.e., the sum over n on the right-hand side 
of Eq. ([50]) . Both convergence and analyticity are crucial for the dcfinitcncss and validity of our proposed approach. 
By analyzing the behavior of the mode sum in the high-frequency limit, we provide here a strong indication this sum 
converges (at least) exponentially in the entire domain r > 2M, and is therefore also analytic in this entire domain. 



1. WKB approximation for the large-o; homogeneous radial functions 



The extended frequency-domain radial functions are give by 

RtmM=Ct m .Rf m M (r>2M), (Dl) 

where the (r-independent) coefficients C lmu) are given in Eq. (f2"Tj) . and {■R^ na) (?")> t( r )l * s a P a i r 01 independent 
solutions to the homogeneous equation 



dri 



V l {r)]R lmui = 0, (D2) 



satisfying suitable boundary conditions at r — > oo and r — > 2M, respectively. In what follows we shall explore the 
asymptotic behavior of the quantities RjL^. at large w, and consequently evaluate the large-w contribution to the 
extended time-domain functions 4>f m {t,r) in Eq. ([55]) . 

Since the potential Vt(r) is bounded (for a given I), in the large-w limit the term in square brackets in Eq. (|D2|) is 
dominated by lo 2 . Using the WKB approximation, we can then write 

Rf m . = [1 - V^r)/^} - 1/A exp Ui J v^a - V^dr^j . (D3) 
The square-root term may be expanded for large uj as 



y/LU 2 -Vi(r) 2£u-Vi(r)/(2u). (D4) 

We shall only consider here the leading-order solution at large w, so we ignore the term Vi{r)/{2uS) in this expression, 
as well as the oc cu~ 2 term in the square brackets in Eq. (|D3[) . We obtain 

RLo = ^ T ' ■ (D5) 
From the form of the ODE (|D2[) it is clear that the Wronskian W is constant, and for the specific pair Rf muJ in Eq. 



(ID5[) it takes the value 

W = 2iu. (D6) 
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2. The source term Zi mLJ and its spectrum 

To calculate the coefficients C^ nuJ we first need to analyze the source term Zi mu (r) and determine its discrete 
spectrum. Inverting the Fourier transform in Eq. (|15|) . we have 

Z lm u{r) = f [1] (r) { r p lm (t, r)e wt dt, (D7) 



where T r denotes the i-period of the radial motion (see App. [A"]) . Throughout this appendix, /I fc l (k — 1,2,3,...) 
denote certain functions of r, independent of w (or n), which are analytic throughout r > 2M (but whose precise form 
would not interest us otherwise). 

Our first goal is to obtain the to spectrum of pi m . For an eccentric geodesic the azimuthal motion may be expressed 

as 

tp p (t) = fl v t+ A<p(t), (D8) 

where £l v is the i-averaged angular frequency dip/dt (see App.|Aj, and A(p(t) is a certain analytic function of t, which 
is periodic with periodicity T r . [Analyticity is directly inherited from that of <p p (t); The T r -periodicity results from 
the fact that for a given orbit dAip/dt is a function of r only] Substituting this in Eq. ([7]) we get 

Pi m (t, r) = Q mg (r 2 t t *)- 1 /3 m («)e-^ t 5[r - r p (i)], (D9) 

where (3 m {t) = e~ jmAv ^^ is an analytic, T r -periodic, function of t. We rewrite this as 

Pim(t,r) = (qfl%(r)p m (t)S[r - r p (i)]) e -' im0 -' 

= S lm (t,r)e- imn ^. (D10) 

The spectrum of pi m is the same as that of Si m , with all frequencies simply shifted by mil v . Now, the dependence of 
Sim on t is only through the T r -periodic functions r p (t) and f3 m {t). Therefore we may write 

S lm (t, r) = J2 Si m n(r)e- inQrt , (Dll) 

n 

and correspondingly 

p lm (t, r) = J2 S lmn (r)e-^ mn - +nn ^\ (D12) 

n 

where tt r — 2tt /T r is the fundamental radial frequency. Thus, the spectrum of pi m (t, r) is the discrete set of frequencies 

lo = m£l(p + nfl r = Lu nm . (D13) 

Finally, we calculate the coefficients Zi mn (r) = Zi muJnm (r) using Eq. (|E)7|) . Pulling the factor /^(r) into the integral 
and using Eqs. ()D10jl and (|D13[) . we get 

Z lmn (r) = q r fl m (r)/3 m (t)S[r - r p (t)]e mn ^dt. (D14) 



3. Calculating the coefficients Cc^ 

We turn now to calculate the coefficients C^ nn = C lmM . These are given by 

= w -r dr> (D15) 

where Rf mn = R^n^ ■ Substituting from Eq. (|D14| for Z\ mn and absorbing the factor l//(r) in //^J (to form another 
analytic function, / ; ^J) we get 

Ctn = iW- 1 f r fl m (r)Rf mn (r)/3 m (t)S[r ~ r p (t)]e w dtdr . (D16) 
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The two integrals in the last equation are allowed to be interchanged (as one may verify, for example, by explicitly 
carrying out first the t- integral and only then the r- integral). One obtains 

Ctn = VW- 1 r fl 4 2K(t)]RLnK(m m (t) e mn ^dt. (D17) 
Jo 

Next we substitute for W from Eq. (|D6[) . and also for using the large-w asymptotic form given in Eq. (|D5|) . 

Denoting r»(i) = r„[r = r p (t)], we thus obtain 

Ctn = / Tr /5r P (i)]fc(i)^ m " mr * (() e l " firt ^ 
Jo 

Jo 

The analytic factor e^™^ 1 "* W (as well as —i/2) may be absorbed in /,[*] — which will in turn become a new analytic 
function, denoted f,J . Defining 

i±=i±r* (D18) 
(namely the two Eddington-Finkelstein null coordinates), the last equation becomes 

Cf mn = q^l r f\TKm m {ty na ^dt, (D19) 
Jo 

where t v ± (t) is t±(r = r p (t))=t± r$(t). 

In the next step we wish to transform the integration variable from t to £§.(t). To this end we briefly discuss the 
properties of this transformation (and its inverse), particularly in terms of analyticity and periodicity. Obviously i±(t) 
is analytic. Also, since the orbit is timelike, i±(i) is monotonically increasing, and dt±/dt nowhere vanishes. This 
implies that the inverse function t(t±) is well-defined and analytic. Note also that all functions of r p a re periodic in 
t±, with the same period T r . The same applies to (3 m (t). Therefore Eq. (|D19jl may be expressed as [3l[ 

GTmn = /^M^IM^)] J^^. (D20) 

The last expression is nothing but (qTrto^) times the n-th Fourier coefficient of the periodic function 

Htnitl) = flTKitl)] PmWl)] A. (D21) 

Since all three factors on the right-hand side are analytic functions of t\, so is H?L Therefore its Fourier 

coefficients must decay (at least) exponentially with |n|, and the same applies to the coefficients C^ mn . 



4. Reconstructing the time-domain extended homogeneous solutions 

We turn now to explore the high-w contribution to the time-domain extended homogeneous solutions 

tLfrr) =E^>K i<W =E C ^n<>) e "' Wmt - (D22) 

n n 

Using the large-u; asymptotic expression. (|D5j) one finds for the large-u; contribution 

n n 
= g-im^t^^g-mfi,.*^ ( D23 ) 

n 

The sum in the last expression has the form of a Fourier series in time with coefficients C^ nn which, as we already 
established, decay at least exponentially. Therefore, this sum converges to an analytic function of t^. It then follows 
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that the entire right-most expression in Eq. (|D23p is analytic in t T . Since t T , in turn, is an analytic function of r and 
t, then so must be this expression. We therefore conclude that, at least within the leading-order approximation, the 
large-ui contribution to (f>i m (t, r) is analytic in both r and t throughout r > 2M. Note also that the sum over n in Eq. 
(|D23[) is guaranteed to converge uniformly for all t and r > 2M. The same applies to the corresponding Fourier sums 
for the r and t derivatives of (j> lm . 

In the above discussion we considered only the leading-order term in the X/u WKB expansion. One naturally 
expects that the contribution from higher-order terms in this expansion will converge even faster, and hence will not 
interfere with the analyticity of r). Also, such higher-order contributions are not expected to affect the uniform 

convergence of the sum over extended n modes. 

We regard the results of the above leading-order calculation as a strong indication that the extended time-domain 
solutions, as they are defined in Eq. (f30|) . are indeed analytic everywhere outside the black hole. This is further 
supported by the numerical results presented in Sec. IIVBI 
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